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Abstract 

The process of downfall of initially homogeneous gas onto a solid ball due to the ball's gravity 
(relevant in astrophysical situations) is studied with a combination of analytic and numerical 
methods. The initial explicit solution soon becomes discontinuous and gives rise to a shock wave. 
Afterwards, there is a crossover between two intermediate asymptotic similarity regimes, where 
the shock wave propagates outwards according to two self-similar laws, initially accelerating and 
eventually decelerating and vanishing, leading to a static state. The numerical study allows one to 
investigate in detail this dynamical problem and its time evolution, verifying and complementing 
the analytic results on the initial solution, intermediate self-similar laws and static long-term 
solution. 



PACS: 47.40.-X, 47.70.Nd, 98.35.Mp 

1 Introduction 



The problem of unsteady motion of compressible gas giving rise to the formation of shock 
waves has been the subject of continuing attention. In particular, the polytropic gas flow in 
one dimension is amenable to powerful mathematical methods [1, 2, 3]. The introduction 
of gravitational forces complicates the problem but widens its range of applications. 

We study in this paper a dynamical problem that is very simple to formulate: given a 
solid body with a spherically symmetric mass distribution (a ball), placed in a homogeneous 
gas, analyze the subsequent evolution of this gas under the body's gravity. This problem 
has applications in several astrophysical situations involving gas accretion. For example, 
it could be a simplified model for the formation of planetary atmospheres: as a result 
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of gas accretion onto the solid (rock) surface of a previously formed spherical core, an 
atmosphere may form. Alternatively, it could represent the downfall of gas onto a neutron 
star. The problem leads to solving the partial differential equations of fluid dynamics in one 
dimension (the radial distance). These equations are nonlinear and no general method of 
solution is available. However, given the simplicity of the initial and boundary conditions, 
many results can be obtained by purely analytic means (as we shall see). The numerical 
integration of the partial differential equations provides an alternative method of study 
which is used here to validate and complement the analytic results. 

Certain aspects of this fluid dynamical problem have been treated before in the as- 
trophysical literature [4, 5]. However, a complete treatment employing the full power of 
nonlinear fluid dynamics methods has not yet been attempted. The thorough analytic and 
numerical treatment provided here allows one to attain a picture of the whole process and 
to connect with other problems in fluid dynamics, especially, with gravitation. In particu- 
lar, the exact description of the formation of the shock wave may illustrate general features 
present in other processes. 

This paper is organized as follows. In section 2, we introduce the relevant magnitudes 
of the problem to get an intuitive idea of the physics involved. Afterwards, we introduce 
the fluid equations to be used. In section 3, an exact solution is found for the initial 
non-trivial dynamics, which occurs near the ball's surface. In the following sections we 
obtain two similarity solutions valid for larger t, the first still confined within a small 
height (section 4), while the second is vahd for large radius (section 5). In section 6 we 
integrate numerically the partial differential equations of the constant gravity and varying 
gravity case and compare the numerical results with the analytical ones. Finally, in section 
7, we consider the long-time asymptotic static state and explore the transition to it with 
the numerical solution. In section 8 we summarize and discuss the results. 

2 Relevant magnitudes and fluid equations 

The initial condition, namely, the (spherically symmetric) solid ball in the homogeneous 
gas, is characterized by four parameters: the radius R and mass M of the ball, and the 
pressure Pq and density Pq of the gas (we assume that the gas is perfect, inviscid, and 
polytropic). In addition, we have the constant of gravity G. From these five dimensional 
characteristic parameters, we can form two independent dimensionless numbers: the ratio 
of densities (4/3)7ri?^(po/^) and the ratio RPq/{pqGM). The quantity Pq/Po ~ Cq (cq 
being the sound speed) is approximately the gas thermal energy per unit of mass. On 
the other hand, GM/R is the gas potential energy per unit of mass on the ball. Their 
ratio measures the relative strength of gravity in our problem. To have any significant gas 
downfall, we must demand that the ball's gravity dominates over the gas thermal energy, 
that is, Rc^/{GM) < 1. 

Defining the scale at which the gas thermal energy is similar to its potential energy as 
TZ — GM/cl, we demand that 
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The last length is the radius of a volume of gas such that its mass is similar to M and, 
therefore, we can neglect the self-gravity of the gas. 

Since the initial and boundary conditions are spherically symmetric, we will assume 
that the whole process is spherically symmetric and thus one-dimensional. We further 
consider an adiabatic evolution of the gas or, more generally, a polytropic equation of state, 
P oc (7 > 1): for a perfect gas, with 7 = Cp/Cy (the constant ratio of specific heats), 
a polytropic process is adiabatic, but the notion of polytropic process is more general. 
For numerical calculations, we will use 7 = 7/5 (the adiabatic index of a perfect diatomic 
gas). Then, we have the continuity equation, the Euler equation and the thermodynamic 
equation: 

dp , d{pv) , 2{fyv) _ 

dt^ dx ^x + R ~ 

dv dv , , 1 dP , , 

dt dx p dx'' 

d P d P ^ 

where x is the radial distance from the ball surface and 51(0;) = GM/ [x+RY = qR^ / {x+R)'^ 
is the gravity acceleration {g is the gravity at the ball surface). The initial conditions are: 
p(0, x) = po, P{^,x) — Pq and v{{),x) = (x > 0). The boundary condition at the solid 
surface is v{t,Qi) = 0. We remark that the behaviour of the solutions of these equations 
for large x is only determined by the initial conditions, which in some sense replace a 
boundary condition "at infinity" . However, in the numerical calculations we shall need a 
real boundary condition at some large value of x. 

The thermodynamic equation (4) (which expresses conservation of entropy or some 
more general thermodynamic quantity) has the trivial solution P/p'' — -Pq/Po- Introducing 
the sound velocity c (c^ = dP/ dp) , we can write the preceding solution as 

c'ip)/cl^ip/poy-' (cl^^Po/po). (5) 
Therefore, we need to consider only the two remaining equations: 

However, in the presence of a shock discontinuity, the thermodynamic equation (4) breaks 
down (as is most intuitive when this equation can be interpreted as entropy conservation). 
Therefore, the solution (5) does not hold, as we shall see in the numerical solution. 

The fluid equations can also be derived as conservation equations, namely, as equations 
for conservation of mass, momentum and energy. In our case, we need to take into account 
that X is the radial distance and write the divergences in spherical coordinates; in particular, 
the divergence of the stress tensor T^^ — pv^v^ + g^^P {g^^ is the inverse metric). The fluid 
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equations in conservative form become: 



dp d[{x + RYpv] ^ 

dt {x + Rfdx ' ^ ' 

djpv) d[ix + R)^pv^] dP _ 

dt ^ {x + R)^dx ^ dx ~ ^ ' 

d{pe + pvy2) , d[{x + Rf{pe + pv^/2 + P)v] _ _ , . 

dt ^ (x + Rydx ~ 

They are equivalent to Eqs. (2,3,4), provided that the pressure is related with the internal 
energy (per unit mass) e by P = (7 — l)pe and that the solutions are smooth. In the 
presence of a shock discontinuity, the conservation equations still hold, which is their main 
advantage. Their disadvantage is that they are more complicated. 

Defining the vector U = (p, pv, pe+pv^/2), we can express the equations in the following 
divergence form: 

with F(U) = {pv, pv'^ + P, {pe + pv^/2 + P)v) and G(U) = {-2pv/{x + P), -2pv'^/{x + 
P) —g{x)p, —2{pe + pv^/2 + P)/{x + R) —g{x)pv) the source term. This form will be useful 
in the next section, and for numerical calculations as well. 



3 Early stage: exact solution 

Initially, the gas will start falling with acceleration g{x), that is, with a negative velocity 
increasing as v ~ —g{,x) t. However, it will be stopped at the ball's surface, where the 
density and, therefore, the pressure must increase. Consequently, a wave transmitting the 
boundary condition f (t, 0) = will propagate outwards. It is therefore crucial to determine 
the law of propagation of waves in the present conditions. This is called, in mathematical 
terms, the analysis of characteristics. Once this is done, and as long as the dynamics 
consists of the propagation of a simple wave, one can obtain an exact solution. We do 
not have here the opportunity to review the theory of simple waves, so we provide some 
indications and refer the reader to general treatises on the theory of one-dimensional gas 
flow for a comprehensive treatment [1, 2, 3]. 

For the moment, we confine ourselves to a spherical shell over the surface of height 
much smaller than P, where we can consider constant g — GM/R?. Hence, the problem 
becomes that of the fall of gas on a flat surface (the ground). So we must use the one- 
dimensional form of the vector divergence and, therefore, we must neglect the last term on 
the left-hand side of the continuity equation (6), writing it as 

The gas still unperturbed by the wave merely falls with velocity v — —g t. In a reference 
system falling with it, it is at rest and, therefore, the speed of sound is cq. Hence, in the 
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original reference system the wave front is located at Xf = cqI — gt^ /'i. Then, for x > Xf, 
the solution is trivial, and we only have to find the solution ior x < Xf. In the falling frame, 
the acceleration of gravity vanishes and the set of two equations (12) and (7) reduces to a 
"gas tube problem" that can be solved by introduction of the Riemann invariants [1, 2, 3]. 

The Riemann invariants are best introduced by writing the equations in the divergence 
form (11). The replacement of spherical geometry with cartesian geometry reduces the 
source term to G(U) = (0, —gpi —gpv) and, furthermore, the change to the falling frame 
removes it. Then we can write the divergence form (11) as a quasilinear partial differential 
equation: 

where A{\J) is the Jacobian matrix of derivatives of F(U) with respect to U. The solution 
of this equation is carried out by diagonalizing A: the eigenvalues A of A (which are real 
because the equations are hyperbolic) define the characteristic curves 

% - (14) 

and the left eigenvectors of A (the eigenvectors of its transpose) define the Riemann in- 
variants J, that is, functions that satisfy 

for some eigenvalue A. It is easy to prove that J is constant along the characteristic 
associated with the corresponding A. The characteristics are the directions of propagation 
of disturbances (waves). Since we have three eigenvalues, namely, v and v ± c, we have 
that the first family of characteristics are the streamlines and the other two families are 
forward or backward characteristics, according to the sign. 

Let us denote a;', f c' the coordinate and variables in the falling frame. Initially we have 
a constant state (with v' and P constant), which is preserved for x' > CqI. Therefore, the 
backward characteristics crossing the forward characteristic x' = c^t transmit a constant 
value of the Riemann invariant J_, so the solution is a forward simple wave. Given that 

2c' 2co 

J- — V — 
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the sound velocity is related with the gas velocity by c' = Cq + (7 — l)v' /2. Furthermore, 
the fact that the solution is a forward simple wave implies that the forward characteristics 
are straight lines (see Fig. 1). Hence, the wave's propagation law v' — F[x' — {v' + c')t] 
is an implicit equation for v', assuming that we can determine the function F. This is 
done using the boundary condition v{t,0) = 0. The implicit equation is algebraic and its 
solution is straightforward. We then obtain: 
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Figure 1: Line described by the ground and forward characteristics in the falhng frame 
(in units such that Cq = g = 1 and for 7 = 7/5). The dashed hne {dx'/dv')t = is the 
envelope of the characteristic hnes. At the cusp point (5/6,5/6), the characteristics cross 
and the solution becomes multivalued. 

For a simple wave, the density is a definite function of the velocity [2]. We can calculate 
it from Eq. (5) to be 



These expressions for p and v constitute the solution of Eqs. (7) and (12) with the given 

boundary conditions. 

Wc note a peculiarity of the solution found above: the wave front's coordinate Xf — 
Cot — begins increasing but, for t > co/g, turns to decreasing and, eventually, returns 
to the origin (a free fall with initial velocity cq). On the other hand, a shock wave singularity 
occurs before the turning point: One can calculate {dx/dv)t, and its null locus defines a 
line in the xt-planc, namely, x = [2co + (7 — l)gt]'^/ (87(7). Initially, the corresponding x is 
larger than Xf, so it is unphysical, but, for t = 2co/[(7 + l)fi']) it meets the wave front; that 
is, the falhng gas overtakes the gas just below it and a shock wave arises. The relevant 
geometry is represented in Fig. 1, in the falling reference frame. From the moment of 
formation of the shock onwards, the solution (16,17) ceases to be valid. Furthermore, the 
analysis of the propagation of a shock wave cannot be done in terms of simple waves. 



When the initial conditions are sufficiently simple, the solution of a one-dimensional gas 
flow problem may adopt a self-similar structure and be greatly simplified by writing the 
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4 Similarity solution with constant gravity 
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equations in the appropriate variables [2, 6], in which the partial differential equations 
become ordinary differential equations (ODE). 

In our case, a self-similar solution is not possible, since we have too many parameters 
in the initial conditions. Nevertheless, it is commonly observed that nonlinear equations 
that do not have similarity solutions at the outset develop them in an intermediate asymp- 
totic regime, that is, a regime between two very different scales where both scales can be 
neglected [7]. We shall see that we can find a similarity solution of this kind. 

In the equations with constant gravity, namely, (3), (4) and (12), R and GM only 
appear in the combination g = GM/ furthermore, the only parameter in the initial and 
boundary conditions is po- Hence, we have one less parameter. However, we can still form 
the two independent dimensionless variables gt/cQ and gx/cl, so we will have to remove 
another parameter to have a similarity solution. 

If we assume that the initial gas temperature is very low, then Cq and the only 
independent dimensionless variable is ,^ = x/{gt'^) (we measure the distance x from the 
ground). This corresponds to the intermediate asymptotics c^/ g ^ x ^ R; equivalently, 
in terms of the dimensionless height variable x* — gx/c^, it corresponds to 1 <^ x* <^ R* — 
n/R (recall that n-> R). 

To take into account the presence of shock waves and the consequent dissipation, we 
consider the full set of equations (3), (4) and (12). We can express them in dimensionless 
form by introducing new variables: 



2 

X „ X 



^ = P = Por, P = — Pop- (18) 



Some straightforward algebra then yields. 



C[u'+{u-2)-]+u = 0, (19) 



r ■ 



{u-2)^u' + -u = -r'^ -r-\2p + ^p'), (20) 
{u-2)^[\n{p/r^)]' + 2{u-l) = 0. (21) 



Upon making the change of variables 

C ~ C ~ 2' 

u = (22) 

(corresponding to changing to the falling reference frame), the term in the second 
equation vanishes and the system of ordinary differential equations simplifies to 

i[u' + {u-2)^]+u = 0, (23) 

(u-2)iu' + -u = -r~\2p + ip'), (24) 
{u-2)i[\n{p/r'^)]' + 2{u-l) = 0. (25) 
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Similarity equations of this type have also been studied in Ref . [4] . However, the authors 
do not make the change to the falling reference frame, and restrict themselves to solving 
the equations by a series expansion near the ground. Instead, we follow here the general 
methods exposed in Ref. [6]. 

In terms of the variable t = In^, the previous equations imply a system of autonomous 
nonlinear ODE, namely. 



du _ P{r) [2 + 7 u{r)] - r{r) u{r) (2 - 3 u{r) + u{rf) 



dr -iPir) + r(r) [u{t) - 2] 



2 



(26) 



^ -2p(r)+r(r) [u{t) - 2] ujr) 

dr (-7P(t) + r(T) [u{t) - 2f) [u{t) - 2] ' 

dp ^ 2 7p(t) - rjr) (4 - (6 + 7) u(r) + 2 u(r)') 

dr ^^^^ -^p(^r)+r{T)[u{T)-2f ' ^ ^ 

Given the homogeneity properties of these equations with respect to p and r, it proves 
convenient to introduce the variable 9 = p/r (the dimensionless form of the temperature 
for the perfect gas with pressure P and density p). Thus, 

d0 _^ 29 (1 + 7(^-2)) -(.2-2) (A-i5 + ^)u + 2e) 

du (m-2) (^ (2 + 7m) -M (2-3m + m2)) ^ ^ 

The solution of this equation provides the relation between the "velocity" u and the "tem- 
perature" 9, 9{u). Substituting it back into Eq. (26), we have an ODE for 'u(t), which is 
immediately solved by a quadrature. Analogously, one solves for r(T). 

We must now translate the initial and boundary conditions of the original partial dif- 
ferential equations (PDE) into an initial condition for the ODE (29): at t = 0, p(0, x) = po, 
P(0, x) = Po = (co 0), and v{0, x) = 0, implying that r = 1 and u = for ^ = C, = 00; 
at a; = 0, f = 0, which implies, since v = gt{^u — 1), that -u = 2 for ,^ = and ^ = 1/2. 
The problem seems overdetermined. However, recalling the exact solution of section 3, 
which gives rise to a discontinuity (shock wave) , we expect that the self-similar solution is 
discontinuous, so that there is no contradiction (in fact, this happens in most self-similar 
solutions [2, 6]). 

When discontinuities arise, the original differential equations break down, but they 
can be expressed as integrals over test functions, which may have discontinuous solutions 
{weak solutions) [1, 2, 3]. Integration of Eq. (11) across a discontinuity yields the Rankine- 
Hugoniot jump conditions F(Ui) = F(U2). These conditions can be written as 

PiVi = P2V2 , 

pivf + Pi = P2V2 + P2 , (r>n\ 
„2 p „,2 p l>5u; 

-?r + hei = — + — + 62, 

2 pi 2 p2 

where subscripts refer to the values on each side of the shock and e is the internal energy; 
for a perfect gas, e = P/[p{l — !)]■ These equations hold in a coordinate system in which 
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the shock surface is at rest. On the other hand, the location of the shock must be at fixed ^, 
say (the subscript s will denote variables pertaining to the shock motion) . Consequently, 
the shock wave velocity is 

dX g , djij n J- r~\ ^ S 

'^^^^ ^^'^ ^ 265^ - 2-, 
that is, Us = 2 (and also Ug = 2). Then, 

ri{ui - 2) = r2(M2 - 2), 

Pi + ri{ui - 2f = P2 + r2{u2 - 2f , . . 

(^1 - 2f ^ 7 Pi _ {U2 - 2f ^ 7 P2 ^ 
2 7 — Iri 2 7 — lr2' 

valid in both the rest and the falling reference frames. Using the falling frame and elimi- 
nating r2/ri between the first and second equations, 

? 6 
^ + -wi - 2 = + U2-2, 



^'-'^ 27 "^'"^ 27 (32) 

(«i - 2)2 + — L-^i = {U2 - 2f + —^^2. 



We have obtained an involutive mapping of the plane (m, d) as the relation between 
the variables at either side of the shock discontinuity. Since we have the condition that 
as a; — > cxo we go to (0, 0) and (in the falling frame) these values hold all the way down 
to the shock surface, we just need the image of the origin under the mapping. This yields 
the point (^4/ (7 + 1), 8 (7 — 1) /(7 + 1)^^ . Therefore, we need the solution of Eq. (29) that 
goes through this point. We can see in the example of Fig. 2 (7 = 7/5) that the strand 
of this solution that departs from (5/3,5/9) ends at the point (2,0), hence satisfying the 
boundary condition at x = 0.^ 

Now, substituting the function just computed back into Eq. (26), we obtain 

-^(m) (2 + 7m) +M (2-3w + 'ii2) ' ^ ' 

from which we can deduce the interval of r between the points (^4/(7 + 1), 8 (7 — 1) / (7 + 1)^ 
and (2,0) by integration. The result for 7 = 7/5 is 0.0880104. Hence, we obtain the quotient 
between the corresponding values of ^. This quotient gives the location of the shock rela- 
tive to the ground: |./(l/2) = exp 0.0880104 = 1.09200 =^ = 0.0919995/2 = 0.0459997, 
so that the coordinate of the shock is Xs — 0.04599975it2. 

Once we have described the self-similar solution, let us analyze in detail its range of 
validity. The self-similar solution is valid for cq — 0, that is, for the dimensionless variables 
t* — gt/co and x* = qx/cq sufficiently large. Moreover, x* <^ i?*, but the upper limit of 
the t* asymptotics is still indefinite. The remaining condition can be found by using the 

^Notice that this imphes that the asymptotics co ^ is of the first kind (the simple case) according to 
the denomination of Barenblatt [7] . 
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0.5 1 1.5 2 

Figure 2: Relevant solution of the ODE (29) for 7 = 7/5, displaying the origin and its 
image at the shock surface, the point (5/3, 5/9). 

dimensionless variables ^ = x*/{t*)'^ — x/{gt'^) and C, = x*/t* = x/{cot), such that the 
former only depends on g and the latter on cq. We have that Cq ^ implies that C ^ 1- 
To realize the consequences of the previous conditions it is convenient to consider the 
logarithmic variables Inx* and Int*, on the one hand, and the linearly related variables 
In^ and In^, on the other hand. The three equations Inx* — Ini?*, InT = and InC = 
define a triangle in the (Inx*, Int*) plane in which the complete intermediate asymptotics 
holds; to be precise, in the interior region far from the borders. Moreover, part of this 
region corresponds to the trivial solution above the shock (the part with ^ > ^s)- 

5 Similarity solution with variable gravity 

For heights of the order of the ball radius or larger, we have to return to the full continuity 
equation (2) and reset g{x) — GM/{x + RY in the Euler equation (3). We have an 
additional parameter, namely, the radius R. If we further assume that i? — 0, we can find 
a similarity solution, that is, in the intermediate asymptotics R <^ x <^ TZ {R* <^ x* <^ 
{R*y). Now, the dimensionless variable ^ becomes ^ = x^ /{GMt^). Hence, we can express 
the continuity and Euler equations as 

r' 

C[Su' +{3u-2)-]+3u = 0, (34) 
r 

{3u - 2) ^u' + u'^ - u = -r-\2p + 3^p'). (35) 

The energy equation becomes 

{3u - 2) e[ln(p/r^)]' + 2{u - 1) = 0. (36) 

These three equations look similar to the ones corresponding to constant gravity, but 
they are more difficult to analyze:^ the change of variables that transformed the latter into 

^An extensive study of self-similar spherical accretion with an initial power-law density distribution is 
provided by Ref. [5]. 



10 



u 



-8 tir 



-0.5 
-1 

-1.5- 
-2 

-2.5 




12 3 4 ^ 

Figure 3: Similarity solution {u{^),r{^)} for free fall with varying gravity. 



an autonomous system is not available. In fact, the free-fall (pressureless) problem is now 
nontrivial, but it can be solved (see Appendix). The solution, in parametric form, reads 

^ 8 cos^ a 

^ 2a + sin(2a)' ^ ' 

, sin(2Q;), sin a 
u = + (38) 

r = , (39) 

Po 9 cos a — cos(3 a) + 12 a sin a 

where ol e (0,7r/2). The corresponding graphs for u and r are shown in Fig. 3. 

As regards the complete similarity solution, it must match the free fall solution outside 
a spherical shock wave with the inner solution with pressure. In turn, this solution should 
be matched with the previous solution near the ground. These matchings arc difficult to 
obtain analytically, so we will rather rely on numerical solutions (next section). At any 
rate, we deduce that, as the shock moves from <S i? to Xg ^ i?, its velocity experiences a 
crossover from = 1x^1 1 = lisgt to = 2/3 {x^jt) ^2/?, {GMQ^/^ t'^/^ (with different 
^s), changing from accelerating to decelerating. In addition, the shock must disappear 
when t ^ GM/cq and Xg^lZ (see section 7). 

Of course, a rough analysis of the range of validity of this similarity solution is also 
possible. It is valid for i? ^ 0, that is, for the dimensionless variables x — x/R and 
i — \/GM / E?/'^ t sufficiently large. It is convenient to use the dimensionless variables 
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^ = x^/i'^ and ( — x/i = ^ R/GM x/t, such that the former does not depend on R. We 

have that — > imphes that C, <^1. Considering the logarithmic variables Inx and Ini, 
on the one hand, and the linearly related variables In^ and In^, on the other hand, the 
three equations ln,x = 0, Int = 3/2 In(7^/i?) and In^ = define the triangle wherein the 
complete intermediate asymptotics holds (far from the borders) . 



6 Numerical solution 

In this section we will numerically validate and complement the main results obtained 
above. As was mentioned in section 2, our aim is to integrate the equation 

for the cases of constant and varying gravity. Relying on the operator splitting approach, 
we integrate first the conservative part of the equation using the Lax numerical scheme 
and then integrate the source term: 



Un+i ^ + _ At F(U-V,)-F(U-_,) ^ ^ ^^^^ 
■10 2 /ax 



With this we get a discrete approximation of the value of U at Xj = jAx and time t = nAt 
to second order accuracy, U". We will show the results in dimensionless variables x* = 

x/{f) and t* = t/{^) and fix Ax* = 0.005 and At* = We have chosen At so that the 
Courant-Friedrichs-Lewy stability condition of the scheme is satisfied [8] . This condition 
esentially states that we are required to choose a time step smaller than the smallest 
characteristic physical time in the problem. At* < , . , which is the time for the velocity 
to lead to a flow over a distance Ax* (during our integration time v* = v/cq < v^^^ = 20). 
We integrate the problem with initial condition p{x,0) = 1, P{x,0) = '^CqPq, v{x,0) = 0, 
adiabatic index 7 = 7/5 = 1.4 and boundary condition at the solid surface v{t, 0) = 0. 

In the case of constant gravity, we can rewrite the equations of gas downfall on a 
flat surface (7) and (12) in the divergence form (40) using G(U) = (0, —gp, —gpv) (as 
was mentioned in section 3). We show in Fig. 4 (left) the density evolution in (decimal) 
logarithmic scale for the case of constant gravity. For the initial stages up to t* = fa 
0.83, before the shock wave develops, the density and velocity values are coincident with 
those given by the analytical solutions Eqs. (16) and (17) (section 3). Then the shock 
formation induces a discontinuity in the density, pressure and velocity distribution and 
the analytic solution is no longer valid. For a polytropic gas and based on the Rankine- 
Hugoniot conservation conditions accross the shock, one can see that the ratios V2/V1 and 
-P2/-P1 can be arbitrarily large, whereas (for these strong shocks) when P2/P1 — >■ 00, the 
density ratio P2/P1 tends to the constant limit P2/P1 = (7 + l)/(7 — 1) [2]- In our case 
with g constant, the pressure and velocity variation across the shock continue to grow as 
long as the shock exists, while the discontinuity in the density tends to the limiting value 
P2/P1 = 6 (7 = 1.4), as is shown in the lateral view of the density in Fig. 4 (right). 
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Following the asymptotic similarity analysis we expect the shock position Xg to move 
as Xg ~ ^s9t^, for the case of constant gravity (section 4), or equivalently x^ — In 
Fig. 5 (left) we plot all the points x*, in the area affected by the wave {x*\p{x* ,t*) > po) 
rescaled by (t*)^. The shock front is at the upper limit of this area. The value ^ = x*/{t*Y 
(evaluated at the shock) evolves to a constant value = 0.046 and therefore, beyond a 
transient time, the shock moves with constant acceleration. This is in agreement with the 
analytic prediction. 

Next we study the time evolution of this process with spherical symmetry for the case 
of varying gravity, with G(U) = {—2pv/ {x + R), —2pv'^/ {x + R) — g{x)p, — 2(pe + pv'^/2 + 
P)/{x + R) — g{x)pv), g{x) = gR^/{x + i?)^ and g the gravity at the surface of the ball. 

For our study case we choose R* = R/ {-^) = 10. We restrict ourselves to this relatively 
small value because the magnitude of the discontinuity at the shock and the velocity of 
the downfaUing gas increase with R*. So the spatial and temporal grid sizes must be fixed 
sufficiently small to be able to handle these sharp discontinuities and satisfy the Courant- 
Lewy condition for maximal velocities in the grid. This in turn imposes a limitation on 
the maximal R* to be considered. 

We remark that, of course, the numerical integration has to be restricted to a finite 
interval of x*, namely [0, x'^^^J, and therefore we need a boundary condition at x'^^^. This 
is not a problem in the constant g case, because what happens for large x* is exactly known: 
the gas is in free fall with acceleration g and p = po (this is a valid boundary condition 
as long as the shock front has not reached a;^^^^). In contrast, in the varying g{x) case, 
we impose as boundary condition the initial values p{t,x*^^^) = po and P{t,x^^^) = Pq — 
7C0P0, which are similar to the real values if » (i?*)^ and as long as the shock front 
is sufficiently far from x'^^^. For our calculations we use x^^^ = 450. 

In Fig. 5 (right) we plot x* /t* vs. t* in (natural) logarithmic scale. The lower area is the 
area between the shock front and the surface. As one can see in the Figure, we can fit the 
motion of the shock front as x*/t* = (t*)"^/^exp (—0.72). We get the corresponding value 

of Cs = (f*)2(ji.*)2 — 0.0011. Over this short range of distances and time, we are already able 
to find the asymptotic behaviour predicted in section 5 in spite of the strong conditions 
needed for the derivation of the similarity solution Cq ^ and R ^ which are hardly 
satisfied in this particular case (for these cases v is not so large as compared to Cq, see Fig. 
6 for a more detailed view of the dynamics of the main physical quantities). 



The dissipation associated with the shock wave makes us expect that the gas must reach a 
stationary state with more entropy than the initial state. The properties of the stationary 
state are easily studied through the hydrostatic equation. 



with P p^ and with either constant or variable gravity g{x) = GMj {x + R^. It is easy 
to obtain the solution in the latter case (the constant g solution can be obtained expanding 



7 Transition to the static state 
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Figure 4: Density distribution and time evolution (for constant g) in (decimal) logarithmic 
scale. Three-dimensional view (left) and lateral view (right) . The density ratio across the 
shock tends to P2/P1 = (7 + l)/(7 — 1) = 6 as expected for a polytropic gas. 



this solution in x/R or directly): 

7-1 



= 1 



7-1 pbQbR 



7 



R 



x + R 



(42) 



where the subscript b denotes values at the ball's surface. It is more convenient to take 
the reference at infinity, where the density and pressure keep their initial values, so that 




Ipo GM 



7 Po x + R 



1 + (7 - 1) 



7^ 



x + R 



(43) 



The transition to the static state is best studied numerically. Next we show, for the 
case of varying gravity with R* — 7, two snapshots of the spatial distribution of the main 
physical quantities {v* — v/cq, the speed of sound in the media c* — c/cq — P*^/ p with 
P* — P/c^i and p/ Po) as obtained from the numerical integration explained in section 6, 
compared with the values of the static solution given by equation (43). In Fig. 6 (left), 
we observe that the shock front appears as a discontinuity of the main variables when the 
speed of the downfalhng gas is greater than the speed of sound (supersonic regime). The 
shock discontinuity decreases in magnitude once the gas speed v* is no longer supersonic. 
In Fig. 6 (right) we can see how the numerical solution tends to the static solution. As the 
shock moves away the gas behind slows down and finally becomes still. 
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Figure 5: Self-similar evolution. We plot, for all the points x* between the surface and the 
shock front, the time evolution of: (left) ^ = x* j (t*)^ for the case of constant gravity and 
(right) [x* /t*) in logarithmic scale for the case of varying gravity with R* = 10. (Left) If 
gravity does not change with height, the value of ^ at the shock front (upper hne) tends to 
a constant value — 0.046. (Right) When gravity changes with height, at the shock front 
(upper line) we can fit x*Jt* = (t*)-^/^ exp (-0.72) and therefore = ^^.y2f^,y2 = 0.0011. 
Both results are in agreement with the analytic predictions of self-similar solutions. 



8 Discussion 

Summarizing, we have divided the space dimension (the radial distance x) into a near zone 
(near the ball's surface, where x <^ K) and a far zone {x ^ R). In the near zone, we 
have obtained an exact solution for short time and a self-similar regime for longer time 
and radius. In the far zone, we have a self-similar regime for relatively long time and 
radius. The parameter controlling the extent of both intermediate asymptotic regimes is 
R* = TZ/R. We have also obtained an exact static state for very long times. Even though 
the regions of validity of the analytic solutions occupy a small part of the plane {x,t), the 
numerical study provides an interpolation between them, in addition to verifying those 
solutions. Therefore, the combined solutions provide a good intuitive understanding of the 
entire dynamics. The crucial feature is the formation and propagation of a shock wave 
until its dissipation. 

We have integrated numerically the equations of this dynamical system with a robust 
scheme that, respecting the flux conservative structure of the equations, is able to reproduce 
the dynamics in the neighbourhood of the shock front. We observe that as soon as the 
gas velocity is above the speed of sound in the media a sharp discontinuity appears in all 
the relevant variables. In the case of constant gravity the magnitude of this discontinuity 
increases up to a constant value whereas the pressure and velocity increase indefinitely. In 
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Figure 6: Simulations for R* = 7: (Left) At t* = 6, in a certain range of x*. the speed of 
the downfalhng gas v* is greater than the speed of sound in the media c* = ^P*j/p and 
a shock front is formed. At this discontinuity the gas is no longer polytropic and P*7/p 
is not given by (p/po)^~^- For comparison we include the asymptotic solution or static 
solution given by equation (43). (Right) At t* — 70 the speed is not supersonic, and the 
gas is polytropic again. As the front moves away from the ball surface, the density and 
pressure tend to the static solution. 



the case of varying gravity, the highest discontinuities are reached in the region Xg <C TZ, 
then the shock front moves away from the ball and the magnitude of the discontinuity 
decreases. It is interesting as well to remark that the shock wave is particularly robust and 
holds for Xg > TZ (see the sharp discontinuity in Fig. 6 (right)). 

Since the propagation of the shock wave mostly takes place in the regimes with neg- 
ligible Pq (if i?* S> 1), it is interesting to consider the one- dimensional fluid dynamics 
equations for pressureless gas, namely, the Burgers equation for the velocity (plus the con- 
tinuity equation, independent of it). The Burgers equation includes viscosity but in the 
limit of vanishing viscosity its only effect is to prevent the velocity from becoming multi- 
valued, giving rise to shocks [9]. Furthermore, the formation of shocks can be determined 
analytically for any initial condition. Therefore, it seems that the intermediate asymptotic 
regimes should be given by the corresponding exact solution of the Burgers equation, which 
is trivial: the gas is in free fall everywhere but the velocity is discontinuous at the ball's 
surface, in order that it be zero on it. In other words, the shock is always at the ball's 
surface and does not propagate. This corresponds to the pressureless dust being deposited 
on the surface, where the density becomes infinite but immediately jumps to its free-fall 
value. Hence, we can appreciate the importance of the dissipation associated to the shock 
wave: whereas in the Burgers equation the kinetic energy of the free-falling dust is simply 
ignored after it inelastically collides with the ground, in the real intermediate asymptotic 
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solutions it serves to heat the gas below the shock wave, producing pressure where there 
was none. This dissipation is connected with the breakdown of the thermodynamic equa- 
tion (4), which, in the case that this equation expresses entropy conservation, indicates 
that entropy is created at the shock. The sharp rise of temperature and pressure behind 
the shock are clearly visible in the numerical simulations. 
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Appendix: Similarity solution for the free fall in the field of a 
point-like mass 

We describe here the "Lagrangian solution" of Eq. (35) with p — 0, that is, 

{3u-2)^u' + -u^ -C^. (44) 

Since the initial condition is given at ,^ — > oo, it is convenient to make the change of 
variable ^ = transforming the equation into 

du ^ u^-u + i 
di i3u-2)i' 

The point (0,0) is a singular point of the ODE. Its solution is unspecified, unless we 

introduce an additional condition. The initial condition v — imphes that v — ^J^^i'^^'^u 

goes to zero as ^ — > 0. To impose it, it is best to linearize and solve the ODE around (0,0): 
The general solution of 

i = (46) 

IS u = + C \fi, so that the singular point is a node. Then, lim|_^Q i~^^'^u = C. Hence, 
the particular solution of the nonlinear equation (45) in which we are interested is the only 
one that satisfies u'{0) = — 1. 

In the Lagrangian picture, the velocity is given by energy conservation as 



v^-^2iE+^). (47) 
The initial condition that the particle is at rest imphes that E — —GM/xq. Hence, 



u^-v^ -. 2^-^(1-^), (48) 

X y xq 
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and to have a self-similar form we must express x/xq as a function of ^. In order to do it, 
we must solve the equation of motion (47): 



_ r = r^^" (49) 



where we have introduced o" = s/xq to write it in a self-similar form. It is clear now that 
= F{x/xq), which can be inverted to obtain x/xq as a function of ^. The substitution 
of this function in Eq. (48) yields the solution of the ODE that satisfies the appropriate 
initial condition. 

One can calculate the integral in Eq. (49) by the change of variables a = cos^ (p. This 
allows one to obtain the solution of the differential equation (44) in parametric form: Let 

cos^a — x/xo; (50) 

then, 

VGMt 1 , sin(2Q;), 

and 

1/2 1 -s/ \ / sin(2Q;), , ^, 

^ = ^cos 3(a)(a+— ^), (52) 

, , , sin(2Q;), . 
u — —cos {a){a-\ — ) sma, (53) 

where a. e (0, 7r/2). It is easy to verify that this parametric solution fulfills the initial 
condition lim^^oow' = — 1. 

In the Lagrangian formulation, the density is given by 



- { ^ 
^ \dx x^ 

where Xq/x^ is a spherical geometry factor. From Eqs. (50), (51) and (52), 



r^L = 8 ^ (54) 

Po 9 cos(a) — cos(3 a) -|- 12 a sin(Qi) ' 

which, together with Eq. (52), constitute the parametric equations of the density. 
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